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The behavior of particles driven through a narrow constriction is investigated in experiment and simulation. 
The system of particles adapts to the confining potentials and the interaction energies by a self-consistent ar- 
rangement of the particles. It results in the formation of layers throughout the channel and of a density gradient 
along the channel. The particles accommodate to the density gradient by reducing the number of layers one by 
one when it is energetically favorable. The position of the layer reduction zone fluctuates with time while the 
particles continuously pass this zone. The flow behavior of the particles is studied in detail. The velocities of 
the particles and their diffusion behavior reflect the influence of the self-organized order of the system. 



I. INTRODUCTION 

Pedestrians in a pedestrian zone [1], ants following a trail 
to food places and many other systems of interacting entities, 
which are moving in opposite directions to each other, show 
a prominent feature, namely the formation of lanes along the 
direction of their motion. This formation of lanes has been 
studied theoretically for colloidal particles in 3 dimensions 
[2, 3, 4] as well as in 2 dimensional systems [5, 6, 7]. These 
examples indicate that flow of particles can have a substantial 
influence on the structure formation of a system of interact- 
ing particles. Experimental studies on such systems have not 
been performed up to date, first hints of a lane formation tran- 
sition could be seen in a 3-dimensional system of oppositely 
charged colloids driven in opposite directions by application 
of an external electric field [8]. Studies of people in panic (for 
example trying to escape from a building) show the influence 
of constrictions on such moving ensembles. 

A system of 2-dimensionally confined moving colloidal 
particles also resembles the classical analogon of a quantum 
point contact in mesoscopic electronics [9, 10] or in metal- 
lic single atom contacts [11, 12, 13]. These contacts exhibit 
transport in electronic channels due to quantization effects. 
Such quantum channels can be seen as similar to the layers in 
the macroscopic transport, since both occur due to the inter- 
action of the particles with the confining potential. A classical 
version of a similar scenario can be built on a liquid helium 
surface, which is loaded with charges. For such a system the 
formation of layers has been reported as well [14]. The change 
of the number of such layers in the vicinity of a constriction 
has been predicted from Langevin dynamics simulations of 
Yukawa particles [15]. 

In biological systems the transport of interacting particles 
through narrow constrictions is of high importance for many 
processes, for example for the size selectivity of transport in 
ion channels [16]. The complexity of such systems allows 
only to make simplified statements on the underlying physics 
governing such phenomena. Experimentally easily accessible 
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model systems can reveal many of the underlying processes. 
In the context of micro-fluidics and "lab-on-a-chip" devices 
one is interested in non-equilibrium transport and mixing phe- 
nomena on the microscopic scale [17]. 

In this paper we present a 2-dimensional system of mov- 
ing, superparamagnetic particles. The interaction energies be- 
tween the particles and therefore the effective temperature of 
the system can be set by application of an external magnetic 
field. The phase behavior of these particles in 2 dimensions 
has been studied extensively [18, 19, 20, 21]. In addition to 
this, it has been shown that confinement of these particles in 
a narrow channel leads to the formation of layers, in order to 
conform to the boundaries set by the hard walls [22, 23]. The 
number and the stability of these layers change as the density 
or the interparticle interactions are varied. In this work we ad- 
dress the question how these layers change when the particles 
are subject to a driven motion along the channel. In order to 
investigate this moving state we first study the properties of 
a static system using Brownian dynamics simulations. Based 
on these results, the moving system is characterized, and the 
results are compared to an experimental system of superpara- 
magnetic particles moving through a lithographically defined 
channel. 



II. EXPERIMENTAL SETUP 

The particles are constricted to a narrow channel connect- 
ing two reservoirs, which are defined on a substrate using 
UV-lithography [24]. Images produced with a scanning elec- 
tron microscope (SEM) of such a channel setup are shown 
in figure 1. Channel geometries of various width and length 
have been produced. The simulation results are compared to a 
channel being 60 \xm wide, 2.7 mm long and having channel 
walls of about 5 \xm in height. The channel is filled with a sus- 
pension of superparamagnetic particles of diameter 4.5 ]im in 
water (Dynabeads). Identical particles have been used previ- 
ously and characterized in [26]. A summary of the properties 
of these colloidal particles is given in table I. 

Gravity confines the particles to the bottom surface of the 
channel due to the density mismatch between the colloids and 
the surrounding water. An external uniform magnetic field 
B = B z is applied perpendicular to the bottom surface. As 
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FIG. 1: SEM images of the full channel geometry connecting two 
reservoirs and an enlargement of the channel entrance/ exit region. 
Also some dried colloidal particles can be seen. During measurement 
the particles outside of the channel are removed so that they don't 
influence the particle transport within the channel. [25] 



diameter a 4.55 =b 0.1 ]im 

mass density /Ww 1-6 g/cm 3 

particle mass m (7.6 =b 0.1) • 10 -14 kg 

saturation magnetization M (5.7 ± 0.4) • 10" 13 Am 2 

effective susceptibility Xeff 7.5 • 10 -11 Am 2 /T 

TABLE I: Particle properties of the Dynabeads used in the experi- 
ment. 

a consequence the colloids form a monolayer in the x-y plane 
with induced parallel dipole moments in z-direction giving 
rise to a purely repulsive pairwise particle interaction. The 
strength of the repulsive force at distance = — r ^ | is 
given by 

V tJ {r i] ) = (^/^)M 2 /r% (1) 

with the magnetic dipole moments M = XefiB of the parti- 
cles. The importance of the pair-interaction can be character- 
ized by the dimensionless interaction strength 

T = /i M 2 (™) 3/2 /(4^£T), (2) 

where n denotes the (overall) particle number density, k b the 
Boltzmann constant, /io = 1.257-10 -6 Vs/Am is the mag- 
netic permeability of free space and T the temperature. For an 
unbounded equilibrated 2D-system which forms a triangular 
lattice, the particle number density can be written in terms of 
the lattice constant a as 



So, T = Vij(a)/ (k B T) is the mean dipolar inter- 

action energy of equation (1) in terms of the thermal energy. 
Accordingly, the applied magnetic field B which is connected 
to the magnetization via M = XefiB plays the role of an ef- 
fective inverse temperature. 



The external magnetic field is the dominant magnetic field 
in this system as it is obvious from the large particle separa- 
tions in the video microscopy snapshot of Fig. 16(a) and the 
mutual induction between the colloids is negligible. Thermal 
and magnetically induced fluctuations of the positions of the 
particles perpendicular to the plane of inclination are less than 
10% of the particle diameter and can be neglected. Tilting of 
the whole channel setup induces transport of the colloids from 
one reservoir into the other due to gravity. An alternative driv- 
ing method would be the application of an in-plane magnetic 
field gradient. 

Before starting experiments the system is set up exactly 
horizontal. The particles are allowed to sediment to the 
bottom surface and arrange in the equilibrium configuration 
within several hours. Before tilting the whole apparatus the 
particles are either all confined in one reservoir (by use of 
laser tweezers) or uniformly distributed along the channel and 
within both reservoirs. In the experiment an inclination of 
c^exp = 0.6° is chosen, where the system is in a gravitationally 
driven non-equilibrium situation, but not yet in the regime of 
plug flow. This inclination results in an average particle drift 
velocity thrift ~ 0.035 (xm/s. A typical snapshot from the 
experiment of the particles moving along the channel is given 
in figure 16(a). 

The particle trajectories are tracked with a video micro- 
scope. The repetition rate of the video microscope setup is 
10 s. All experiments are made at room temperature T w 
295 K. 

In the experiment the number density of the colloids is de- 
fined as the number of colloids divided by the area of the 2D 
channel within the field of view of the video microscope ac- 
cessible to the centers of the colloids. This dimensionless 
parameter T was introduced by Zahn et ah [18], who stud- 
ied experimentally the so-called KTHNY phase transition in 
an unbounded two-dimensional equilibrium system of super- 
paramagnetic particles, to characterize the system state. They 
found that for T < = 52.9 the system behaves like a fluid, 
and for T > T m = 60.9 the system forms a triangular lat- 
tice. For the T- values in between they observed the so-called 
KTHNY or hexatic phase. 

In the experiments described below a magnetic field of 
strength B = 0.24 mT is applied, corresponding to T w 72 
which is in the solid state region of the phase diagram. 



III. SIMULATION DETAILS 

We conduct Brownian dynamics (BD) simulations of a two- 
dimensional microchannel setup in order to investigate the 
flow behavior of the colloidal particles within the channel sys- 
tematically for various parameter values of inclination, over- 
all particle density, and channel width. The equation of mo- 
tion for an individual colloidal particle is given by an over- 
damped Langevin equation. This approach neglects hydro- 
dynamic interactions as well as the short-time momentum re- 
laxation of the particles. Both approximations are fully jus- 
tified in the current experimental context. Typical momen- 
tum relaxation times are on the order of 100 |us and therefore 
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much shorter than the repetition rate of the video microscopy 
setup (10 s) used in the experiment. Thus the colloidal trajec- 
tories Ti(t) = (xi(t), Vi(t)) (i = 1, . . . , N) are approximated 
by the stochastic position Langevin equations with the Stokes 
friction constant £ 



C^f 1 = -V P4 J2 Vij{r iS ) + FT + Ht). (4) 



The right hand side includes the sum of all forces acting on 
each particle, namely the particle interaction, the constant 
driving force along the channel F® xt = mgsm(a)H and the 
random forces Fi(t). The latter describe the collisions of 
the solvent molecules with the ith colloidal particle and in 
the simulation are given by a Wiener process, i.e. by ran- 
dom numbers with zero mean, (F^(£)) = 0, and variance 
(Fia(t)F i(3 (0)) = 2k B T£5(t)5 ij 5 a(3 . The subscripts a and 
(3 denote the Cartesian components. The effective mass m of 
the particles is determined by the density mismatch between 
the particles and the solvent. These position Langevin equa- 
tions are integrated forward in time in a Brownian dynamics 
simulation using a finite time step At and the technique of 
Ermak [27, 28]. 

Particles are confined to the channel by hard walls in y- 
direction and at x = (channel entrance). These walls are re- 
alized both as ideal elastic hard walls and as proposed in [29], 
where a particle crossing the wall is moved back along the 
line perpendicular to the wall until contact. Both realisations 
result in the same flow behavior. Also we performed simula- 
tions with the particles at the wall kept fixed. The channel end 
is realized as an open boundary. To keep the overall number 
density in the channel fixed, every time a particle leaves the 
end of the channel a new particle is inserted at a random po- 
sition (avoiding particle overlaps) within the first 10% of the 
channel, acting as a reservoir. A cutoff of 10a was used along 
with a Verlet next neighbor list [28]. Checks of particle over- 
laps are included in the simulation, but for all ordered systems 
we never found two overlapping particles. 

Starting from a random particle distribution within the 
channel, we first calculate an equilibrium configuration 
(pext _ o) of a closed channel with ideal hard walls. Af- 
terwards we apply to the configuration of uniform density the 
external driving force and allow the system to reorganize for 
10 6 time steps, before we evaluate the configurations. The 
time step At = 7.5 • 10" 5 t# is used, with r B = ^a 2 /k B T 
being the time necessary for a single, free particle in equilib- 
rium to diffuse its own diameter a. We choose £ = Zixr)o, 
with r\ denoting the shear viscosity of the water. The sim- 
ulations are done with 2000 — 4500 particles, for a chan- 
nel geometry of L x = 800a and L y — (9 — 12) a, and 
Xeff,sim = 3- 10 -11 Am 2 /T. Thus external magnetic fields 
B = 0.1 — 1.0 mT and a total particle density of n = 0.4a -2 
correspond to T ^ 21.34 - 2134. 



IV. EQUILIBRIUM PROPERTIES OF THE CHANNELS 

Equilibrated configurations of systems confined to a mi- 
crochannel are used as starting configurations for our analysis 
of the transport behavior. This guarantees that at the begin- 
ning of the transport simulation the particles are uniformly 
distributed over the whole channel. First, we compare some 
results found for the 2D microchannels in equilibrium (the ex- 
ternal driving force is switched off) with the results of Hagh- 
gooie and coworkers [22, 23, 30]. 

During the equilibration process the channel beginning at 
x = and the channel end at x = L x are either closed by ideal 
hard walls, or periodic boundary conditions are applied in in- 
direction. By doing so, we assure that no transport is initiated 
due to the boundary conditions used. The simulation start pa- 
rameters are chosen in such a way that they closely reflect the 
situation of the experiment. In all simulations the area L x • L y 
is defined as the region accessible to the particle centers. This 
is the reason, why in the following simulation snapshots the y- 
positions of the edge particle centers coincide with the channel 
boundary. When comparing the channel widths in the simu- 
lation to the widths of the channel in the experiment, one has 
to add the particle diameter a resulting in L y xp = L y + a, 
e.g. a channel with L y = 10a corresponds to a channel of 
L y xp — 11a = 49.5 ]Ltm for the particles used. The equilibra- 
tion process is usually started from a uniform random particle 
distribution over the whole channel. But to avoid a physical 
instability of the starting configuration the particle separations 
are limited to values greater than 0.7<j. For very dense systems 
this initialization method of course breaks down and we start 
from a hexagonally ordered configuration. 



A. Influence of the Confinement 

The triangular lattice is the high density equilibrium con- 
figuration of an unbounded 2D system. Here, we analyze how 
the confinement modifies the resulting equilibrium configura- 
tions. We submitted simulation runs to determine the equilib- 
rium configuration in dependence of the channel width L y for 
a superparamagnetic system with B = 0.5 mT applied and 
the global particle density n = 0.4a -2 which corresponds to 
T = 533.74 and is deep in the solid phase region. Typical 
snapshots of representative parts of the equilibrium configu- 
rations being obtained are shown in Fig. 2. Shown are the 
regions 300a < x < 600a of a channel with a total length of 
L x = 800a. Notice, that the channel widths are stretched by 
a factor of about 6.67. 

Obviously, whether an ordered or a perturbed configuration 
is formed strongly depends on the channel width L y . For cer- 
tain channel widths it is energetically favorable for the system 
to arrange into what we call layers. Right of configuration 
snapshots of Fig. 2 the equilibrium density profiles transverse 
to the channel walls are plotted. They are calculated by tak- 
ing the average over 2000 equilibrium configurations. For the 
channel widths L y = 7a, 8a, and 10a the peaks of these den- 
sity histograms are well separated and occur at almost regular 
spacing across the channel. These properties are the signature 
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L y = 20a a clearly boundary induced layered structure occurs 
for a system at T = 133.44. This is shown in Fig. 3. 



B. Layer Order Parameter 

The number of layers forming within the channel can be 
identified by an appropriate local order parameter. We there- 
fore divide the channel of width L y into several bins in re- 
direction each containing nt>i n particles and define for differ- 
ent number of layers n\ the so-called layer order parameter 



layer, n\ 



1 

Wbii 



"™ . 27r(n,-l) 



(5) 



which is unity for particles distributed equidistantly in n\ lay- 
ers across the channel width starting at y = 0, and signifi- 
cantly smaller for the non-layering case. 



FIG. 2: Typical simulation snapshots of partitions with length 300cr 
of equilibrated configurations for a dipolar system (B = 0.5 mT, 
T = 533.74) and a selection of channel widths (10a, 9a, 8a, 7a, 
and 6a from top to bottom). The channel widths are stretched by 
a factor of about 6.67. All configurations have the overall particle 
density n — 0.4a -2 . The red curves at the right of each configu- 
ration snapshot show averaged density profiles across the channel. 
For clarity reason, the large magnitude peaks at the walls have been 
truncated at a fixed peak height. 
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FIG. 4: Simulation: Comparison of the local layer oder parameter 
flayer, n z =7 (x) for 7 layers and the local orientational order parame- 
ter ^q{x) along the channel. 



of a well defined layered structure parallel to the walls. For 
the channel widths L y = 6a and L y = 9a the system cannot 
equilibrate into such a single layered structure over the full 
channel and only partial layering is visible in the configuration 
snapshots. Such a confinement induced layering phenomenon 
is in agreement with the results for liquid-dusty plasmas [31] 
and the results of the simulations of Haghgooie [22] . 

The channel widths of 10a, 9a, 8a, 7a, and 6a correspond 
to the widths 6.80R, 6.12R, 5.UR, 4.76R, and 4.0SR in units 
of R = 1.471a, which is the expected separation of layers 
for the unbounded system. Even for wide channels of width 
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FIG. 3: Simulation: Full density profile transverse to the confining 



An exemplary comparison between the results of the local 
layer oder parameter ^iayer,n z =7(^) for 7 layers and the local 
orientational order parameter 



(6) 



which returns a measure of the orientational order based upon 
the distribution of angles 6^ (measured with respect to a fixed 
axis) of the lines joining a particle i with its surrounding 
neighbors, is shown in Fig. 4 for a channel of width L y = 10a 
and T = 533.74 as depicted in Fig. 2. Both order parameters 
have been averaged over 500 equilibrium configurations. The 
equilibrium system consists of 7 layers which is indicated by 
^iayer,n z =7(^) having values close to unity. The small off- 
set results from not fully equidistant peak separation. The 
distance between the central layers is slightly greater than 
the distance between the wall layer and the layer next to the 
wall. The local orientational order parameter ^(x) has val- 
ues greater than 0.6, the signature of a nearly triangular sys- 
tem, but exhibits several dips along the channel length. These 
are connected to the occurrence of defects, i.e. bulk particles 
having 5 or 7 nearest neighbors instead of six and edge parti- 
cles with 3 or 5 nearest neighbors. The nearest neighbors of 
each particles are determined by a Delaunay triangulation. 
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System State Dependency on the Channel Width 



FIG. 5: Simulation: Snapshot of a partition (450a < x < 580cr) of 
the equilibrium defect configuration for the system with L y = 10a 
as shown in Fig. 2. Full circles (•) mark the bulk particles with 6 
nearest neighbors and particle on the wall with 4 nearest neighbors, 
symbol x corresponds to fivefold symmetry (or threefold if on the 
wall), and symbol V to sevenfold symmetry (or fivefold if on the 
wall). 



In Fig. 5 all defects within a partition of the equilibrated 
configuration are marked. For the layered system state the 
defects always occur in pairs (forming a dislocation) and are 
located predominantly close to the walls with quite a regular 
spacing. Due to the purely repulsive nature of the particle pair- 
interaction the edge particles are pressed against the confining 
ideal hard walls as it is obvious from the high peaks of very 
small width at the boundary of Fig. 2. These defects along the 
walls are a consequence of a (slightly) higher line density of 
the edge particles compared to the bulk layers. For example, 
for the system with L y = 10a of Fig. 5 the line density of the 
wall layers is about 6% higher than of the nearest bulk layers. 
Edge layers have only a single neighbor layer whereas bulk 
layers have two. Putting an additional particle into a layer 
results both in stronger interaction within this layer and of 
this layer with its neighboring layers. Thus, it is energetically 
favorable for the system to have defects along the wall instead 
within the bulk, because there the involved energy barrier is 
lower. 

The appearance of dislocations along the wall was also seen 
in [32, 33], where we systematically analyzed the equilibrium 
configurations constricted within a circular hard- wall confine- 
ment for dipolar and screened Coulomb pair interaction as 
function of the particle number. In these systems the parti- 
cles arrange in multiple circles and the defects occur due to 
the bending of the lattice in presence of the curved boundary. 
This is in contrast to the situation here, where the planar walls 
give no need for the lattice to bend. 

So, we can conclude that the layer order parameter is more 
suitable than ipe(x) for the detection of layered structures and 
changes therein, because it is insensitive to defects close to the 
wall. 



C. "Phase Diagram" of the Laterally Confined Dipolar System 
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FIG. 6: Simulation: The layer order parameter as function of the 
channel width. The simulation parameters are: B = 0.25 mT, V = 
133.44, R = 1.471cr, L x = 800cr, and periodic boundaries in re- 
direction. 

The influence of the channel width on the system state is 
analyzed by examining the behavior of the global layer or- 
der parameters ^i aye r, m • The result is shown in Fig. 6 for 
channel widths between 2a and 10a. The global layer order 
parameters as function of the channel widths show for differ- 
ent number of layers n\ distinct response regimes where their 
values are close to one. On top of the graph we also indicated 
the channel width in units of the length scale R. Clearly, the 
change of the number of layers happens with a period of ~ R. 
But for integer multiples of R the system is not in a layered 
configuration, but in the transition between two layered struc- 
tures. This means that the confinement induced optimal layer 
separation is smaller than the separation R expected for the 
unbounded system. 



2R 3R 4R 5R 6R 7R 




I i i " i " i " i " i " i " i " i " I 

23456789 10 11 
L y [a] 



Two independent simulation parameters have a strong in- 
fluence on the state of the dipolar system laterally confined 
between two parallel ideal hard walls. These are the wall 
separation L y and the dimensionless interaction strength T. 
In the following we will compare these dependencies for our 
simulation parameters qualitatively with the results of Hagh- 
gooie [30]. 



FIG. 7: Simulation: The bulk defect concentration as function of the 
channel width for identical simulation parameters as in Fig. 6. 

The above scenario can be confirmed by looking at the bulk 
defect concentration 

Cb — defect sn\ 

defect = jy b \') 
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which is defined as the ratio of the number A^ efect of bulk 
particles with either more or less than six nearest neighbors 
and the total number N b of bulk particles. All particles with 
a distance greater then 0.5a are defined as bulk particles. In 
Fig. 7 C\ eiect is plotted as function of the channel width for 
identical simulation parameters as used above. The concen- 
tration of defects in the bulk shows an oscillatory behavior 
with a period of ~ R. The peak positions indicate the chan- 
nel widths where the system can not equilibrate into a layered 
structure, and the positions of the minima coincide with stable 
layer configurations. This behavior is in good agreement with 
the results of Haghgooie as can be seen from taking slices of 
constant Yh in figure 6 of [30]. 
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Time Evolution of the Defect Configuration 



FIG. 9: Simulation: Density profiles transverse to the walls for L y = 
9a and L y — 10a in dependence of V. Again, the peaks at the walls 
are truncated for better clarity. 
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FIG. 8: Simulation: Time evolution of the defect concentration 
^defect for a channel with n = 0.4a" 2 , N = 4000, and L y = 10a 
kept fixed. The particle interaction strength is modified via the ap- 
plied magnetic field B having values between 0.1 mT and 0.5 mT. 

In Fig. 8 the time evolution of the defect concentration 
C'defect of the bulk particles during an equilibration run is ex- 
plicitly plotted for a selection of T values for a channel of 
width L y = lOcr. All runs are started from a random particle 
distribution. After a time of 10 tb the defect concentration re- 
mains unchanged for all Y values. For 45.0 < T < 80, i.e. for 
the transition region between the liquid and the solid state, the 
equilibration process is slower than for the other values. The 
fluctuations increase near the phase boundary. These effects 
are consistent with the results of Haghgooie [30] obtained for 
an unbounded system. 



System State Dependency on the Interaction Strength 



In Fig. 9 we show density profiles transverse to the confin- 
ing walls for the two channel widths L y = 9a and L y = 10a 
at four values of T. On the left hand side both systems are 
liquid whereas on the right hand side they are both in the 
solid state. These density histograms are obtained by tak- 
ing the average over 3500 configurations in equilibrium. The 
system characteristics are very different depending on the T 



value and the channel width L y . For high Y values, where the 
system is in the solid state, the density profile for the chan- 
nel width L y = 10a is sharply peaked at the positions of the 
seven layers. On decrease of the interaction strength Y these 
peaks broaden and have a Gaussian profile down to a value of 
r « 65. The central peaks show greater broadening than the 
peaks at the wall, i.e. the system melts first in the center of the 
channel. Even for low Y values as Y « 12.01, where the un- 
bounded system would be deep in the liquid state, the particles 
at the wall are still relatively localized in their ^-positions. A 
clear density minimum between the colloids in the edge layer 
and the colloids of the central region can always be identified. 
For the channel width L y = 9a the melting scenario is dif- 
ferent. The peak profile is less pronounced for Y = 533.74 
and there is less order across the channel. A mixture between 
a structure of 6 and of 7 layers is indicated by the positions 
of the peak maxima. The structure of seven layers is favored 
more, because the peaks connected to a structure of 7 lay- 
ers are more pronounced than the remaining peaks related to 
6 layers. Decreasing Y again leads to a broadening of the 
peaks and the structure with six layers becomes more favor- 
able (r = 133.44). The unbounded system would be well in 
the solid state at this value at this interaction strength. For 
T = 85.40 only the peaks related to six layers remain, and for 
T = 12.01 no significant qualitative difference to the situation 
for the channel of width 10a exists. 

These changes of the peak characteristics of the density pro- 
file across the channel of width L y = 9a is also reflected 
in the behavior of the layer order parameters in Fig. 10 for 
n\ — 6 and n\ = 7 layers on variation of the interaction 
strength. flayer, n z =6 exhibits a maximum at about Y = 90, 
and strongly decreases for higher Y values whereas the values 
of ^ layer, m=7 increase to values of about 0.8. 

In Fig. 11 the behavior of the bulk defect concentration 
C'defect anc * the layer order parameter ^layer, m on variation 
of r are summarized for a selection of channel widths. The 
curves are color coded depending on whether the equilib- 
rium configuration has a boundary induced layered structure 
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FIG. 10: Simulation: Comparison of the dependency on the inter- 
action strength F of the global layer order parameter ^i aye r, n t with 
ni — 6 and ni — 7 for the channel width L y = 9. 
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FIG. 1 1 : Simulation: Order parameters in dependency of the dimen- 
sionless interaction strength T. In (a) the bulk defect concentration 
^defect and in (b) the layer order parameter ^i aye r, n t are shown for 
a selection of channel widths. 



(red curves) or not (green curves). The blue curves are con- 
nected to the channel width L y = 6.1<j, where the equilib- 
rium system has a perturbed structure with 5 layers as it can 
be deduced from ^i ay er,n z =5 > 0.9 and C^ efect « 0.25 for 
T > 300. Particles changing between the central and its two 
neighboring layers perturb the 5 layers. 

The red curves for the defect concentration of bulk parti- 
cles C^ efect in Fig. 11(a) show a very similar behavior. All of 
them monotonically decrease up to T « 100 to values < 0.1 
and stay constant thereafter. Only for the small channel width 
L y = 5a which has an equilibrium configuration of 4 layers 
the final defect configuration is about 0.17. Due to the small 



channel width the defects being induced in the layers next to 
the edge layers (because of the higher line concentration of 
the edge particles) are of greater influence. For L y > 6.0a the 
blue and the green curves of C^ efect (T) also show a monotonic 
decay, but of varying magnitude and the C\ eiect becomes con- 
stant at significantly higher T values than for the red curves. 
It is interesting to note, that the curves for L y = 6.0a (green) 
and L y = 6.1a (blue) fall on top of each other for T < 105, 
but significantly diverge for T > 105 where layers form for 
L y = 6.1a but not as strong for L y = 6.0a (cf. Fig. 11(b)). 

In Fig. 11(b) the global layer order parameter ^i a yer, m is 
plotted as function of the interaction strength T. Shown are 
the functional dependencies of ^i a yer, m for the parameter n\ 
which have the maximum value for T > 500. The layer order 
parameters ^i a yer, m increase monotonically to values greater 
than 0.9 for the systems connected to the red curves. The 
transition to the layered structure takes place for T < 100. 
In general, Fig. 11(b) shows that in case of layer formation, 
larger values of L y require larger T values for layering. As 
observed before in Fig. 10 the ^ia y er, m have non-monotonic 
behavior and the transition to the final state takes place for 
T > 110, which is greater than for the layered structures. The 
highest values of ^i a yer, m are less than 0.9. 

Piacente and coworkers [34] studied the structural, dynami- 
cal properties and melting of a quasi-one-dimensional system 
of charged particles, interacting through a screened Coulomb 
potential in equilibrium. This system is related to our situa- 
tion, but a different particle interaction potential is used and 
the particles are confined in ^-direction by a parabolic poten- 
tial. They also find a rich structural phase diagram with differ- 
ent layered structures as function of the screening length ft^ 1 
and the electron density n e of the system. 

A re-entrant phase behavior, i.e. a melting process suc- 
ceeded by a system solidification and subsequent further melt- 
ing, was observed for particle confinement inside of a cir- 
cle [26, 32, 35, 36] or in static ID periodic light fields [37, 
38, 39, 40] both in experiment and simulation. For our pla- 
nar wall confinement we do not find any re-entrant behavior 
as function of the dimensionless interaction strength T (the in- 
verse effective temperature). In Fig. 11(a) the defect concen- 
tration of the bulk decreases monotonically with increasing Y 
and thus gives no hint on a reentrant behavior. This observa- 
tion again is in agreement with the results of [30]. 

For particles inside a disc shaped cavity the increase of ra- 
dial fluctuations is responsible for the re- stabilization of an 
ordered shell structure with increasing temperature. In our 
case the influence of the confining hard walls does not seem to 
have a similar effect on the particle fluctuations in ^-direction 
to give rise to a re-entrance behavior. We conclude that the 
re-entrance phenomenon depends strongly on way of confine- 
ment. It would be interesting to study the influence of the 
curvature of the confinement on the melting scenario system- 
atically. On the other hand, a boundary induced reentrant be- 
havior between different layered structures is observed for in- 
creasing channel width L y (cf. figures 6 and 7). 

Figure 12 illustrates the qualitatively different particle mo- 
bilities due to the confinement according to their ^-position 
for the two channel widths L v = 9a and L v = 10a. Shown is 
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FIG. 12: Simulation: Superimposed particle positions in equilibrium 
for channel segments of length 60a and widths L y = 9a and 10cr 
during the time interval At = 15tb (= 2 • 10 5 BD steps). Shown is 
the situation (a) in the fluid regime (T = 56.38 or B = 0.1625 mT) 
and (b) in the solid state (T = 133.44 or B = 0.25 mT). For the 
images of (b) the corresponding density profiles are given in Fig. 9. 



the overlay of 1000 equilibrium configurations, which corre- 
sponds to a run of length At = 15r#, both for the fluid state 
[Fig. 12(a)] and the solid state [Fig. 12(b)]. We recognize al- 
ready from these superimposed snapshots that the walls affect 
the particle mobilities transverse to the walls for both widths. 
The edge particles have a very low mobility to move away 
from the confining walls they are pressed against. A clear 
depletion zone exists between the edge and its neighboring 
layer. Generally, the spreading of the particle positions in y- 
direction increases with growing distance to the walls. When 
comparing the two widths L y = 9a and L y = 10a we again 
realize that the layered system of width L y = 10a is higher 
ordered with smaller spreading of the particle positions. In 
the fluid state, shown in Fig. 12(a), the particles still move 
predominantly within the layers parallel to the walls. This 
boundary induced layering effect is stronger for L y = 10a 
than for L y = 9a. 

The effect of the type of confinement on the ordering of a 
crystal confined to stripes of finite width was analyzed using 
Monte-Carlo simulations by Ricci and coworkers [41, 42, 43]. 
In their case, the particle pair interaction is given by the in- 
verse power law ex r -12 . They studied the influence of ideal 
planar hard walls and structured walls obtained by fixing the 
wall particles at separations they would have in a bulk system. 
Our findings are in good agreement with their results. 



D. Diffusion Behavior 

For channels, which are small enough, so that the particles 
cannot pass each other the diffusion behavior of the particles 
changes. The sequence of the particles remains unchanged 
and the particles move in a single file (SF). The long-time be- 
havior of the mean- square displacement for infinite long chan- 
nels is predicted to be [44, 45] 

{Ax 2 ) = 2FVi. 



Here F is the single file mobility and t the time. 

Such a behavior is an example of anomalous diffusion 
or non-Fickian diffusive behavior, which is characterized by 
the occurrence of a mean- square displacement of the form 
(Ar 2 ) oc t a , where a ^ 1. The motion is called sub-diffusive 
for the (anomalous) diffusion coefficient < a < 1 and 
super-diffusive for a > 1. The phenomenon of single file dif- 
fusion (SFD) has received a lot of attention in recent publica- 
tions, especially after the experimental observation of Wei et 
al [46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59]. 
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FIG. 13: Simulation: Single file diffusion behavior for a channel 
setup with L y — 0.5a, n — 0.8a -2 , V — 60.39, and periodic 
boundary condition in ^/-direction in absence of any driving field. 
The two symbol types used for the calculated data points refer to two 
different ways of evaluation of the MSD. 

In Fig. 13 we plot the mean square displacement (MSD) 
as a function of the simulation time in a double logarithmic 
graph. The data points are obtained from a simulation run 
of a channel having ideal hard walls, the width L y = 0.5a, 
and periodic boundary condition in x-direction and no driv- 
ing field. Two different algorithms have been used to evaluate 
the MSD. Both, the conventional analysis of the MSD (green 
crosses) and the so-called order-n algorithm (red squares) to 
measure correlations being introduced in [60] give identical 
results. At short times, i.e. at times less than 0.1r#, the MSD 
increases ex t, which is characteristic for the ballistic move- 
ment of the particles. The dashed magenta line in Fig. 13 has 
the slope a = 1 as it is the case for normal diffusive behav- 
ior. Clearly, for t < O.Hb the simulation data points fall 
onto this curve. After the time t#, which can be interpreted 
as the time a particle needs to meet one of its nearest neigh- 
bors and to realize it cannot overtake, the MSD approaches 
the square-root time dependency characteristic for SFD. This 
is indicated in Fig. 13 by the solid blue line, which is a fit 
of the function f(t) = A-t a with the two fit parameters A 
and a to the data points with t > tb> The resulting slope is 
a = 0.5022 ±0.0048, which is in perfect agreement with SFD 
behavior. 

Now, the following questions arise: How does the longitu- 
dinal and transversal particle diffusion behavior depend on the 
channel width L y l How does the transition take place from 
(8) the single-file diffusion behavior for channel widths where 
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particle can not pass each other to the Fickian diffusion be- 
havior of bulk systems? 
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FIG. 14: Simulation: Comparison of the mean-square displacement 
(a) (Ax 2 ) parallel and (b) (Ay 2 ) perpendicular to the confining 
channel walls for the two channel widths L y = 9a and L y = 10a. 
The simulation parameters are: L x = 800cr, n = 0.4a -2 , B = 
0.2 mT, and F = 83.4. The periodic boundary condition is applied 
in x-direction. 

To give a first answer to these questions, we evaluated the 
time-dependency of the mean-square displacement for the two 
channel widths L y = 9a and L y = 10a. The results of the 
MSD (Ax 2 ) parallel to the channel walls and of the MSD 
(Ay 2 ) transversal to the channel walls are shown in the fig- 
ures 14(a) and (b). Obviously, the dependency of the MSD 
on time differs strongly for the two channel widths. The long 
time behavior of the longitudinal MSD (Ax 2 ) is linear propor- 
tional to time t for the width L y = 9a, whereas for L y = 10a 
the long time behavior scales with the exponent a « 0.67. 
For times t < OAtb the longitudinal MSD scales approxi- 
mately with the exponent a = 0.5. For the transversal MSD 
(Ay 2 ) the time dependency is similar (cf. Fig. 14(b)), but for 
the width L y = 10a the transversal MSD has a plateau for in- 
termediate times. This is due to the boundary induced forma- 
tion of seven layers, which is not the case for width L y = 9a. 
Notice, that for either (Ax 2 ) or (Ay 2 ) the absolute values in 
the intermediate and long time limit are smaller than it is the 
case with well defined layers. 

The crossover from single-file diffusion with the exponent 
a = 0.5 to Fickian diffusion with a = 1 in the bulk limit 
has recently been analyzed theoretically by Mon, Percus and 
Bowles [61, 62, 63, 64]. These authors present a phenomeno- 
logical theory in terms of the hopping time Th op , which is de- 
fined as the average time a particle must spend before it can 
"hop" over (pass) its nearest neighbor in longitudinal direc- 
tion. They theoretically show that with increasing transversal 
system size the diffusion constant will increase from zero ac- 
cording to D oc (t^p) -1 / 2 , where Th op is a function of the 
pore radius 3D or the channel width in 2D. They confirmed 
this predicted behavior by MC and MD simulations of hard 



spheres within a pore and hard discs confined to a microchan- 
nel respectively. In general, when particles are allowed to pass 
their neighbor particles, the long time dynamics is given by 
Fickian diffusion. This can be understood by the simple ar- 
gument: After the mean time Th op a particle passes one of 
its nearest neighbors in either direction. Therefore the long 
time diffusion behavior is given by conventional Fickian dif- 
fusion, but for times less than Th op the SFD behavior is ex- 
pected. A similar argument was already used in the context of 
a two chain lattice gas model of Kutner et al. [65]. 



V. TRANSPORT BEHAVIOR OF COLLOIDS IN 
MICROCHANNELS 

Now, we want to address the transport behavior of colloids 
confined to such microchannels as described in the previous 
section. The colloids are driven by the application of an exter- 
nal driving force F and thus form a system in non-equilibrium. 
This driving can be of gravitational origin as in our case, 
or due to the presence of an electrical or magnetic field or 
an osmotic pressure difference between both channel ends. 
To match the experimental situation closely, we will concen- 
trate mainly on colloids with repulsive dipolar pair-interaction 
driven by gravity. First we introduce the effect of dynamical 
rearrangement of the colloids during their flow along the chan- 
nel. We call this effect layer reduction. 



A. Layer Reduction 

A first impression of the particle arrangement under the in- 
fluence of an external driving field give the figures 15 and 16 
which depict typical configuration snapshots from simulation 
and experiment. The particles move along the channel from 
left to right in the positive x-direction. The external magnetic 



10 




k^H^MM^nniHHlHHiHrinH# 

100 200 300 400 500 600 700 800 

x[tr] 



FIG. 15: Simulation: Full channel snapshot for a channel with ideal 
hard walls (T = 533.74, a = 0.2°) after 10 6 BD simulation steps 
having reached a stationary non-equilibrium state. Note, that the 
scaling on the y-axis is stretched by a factor of 20 compared to the 
x-axis scaling. 

field strength B, which is responsible for the strength of the 
pair-interaction, and the overall particle number density n are 
chosen in such a way, that the confined equilibrium system is 
hexagonally ordered. This is true also for the unbounded sys- 
tem under identical conditions. Figure 15 is a representative 
snapshot taken in the simulation of the full channel having the 



length L x = 800<r = 3.6 mm. The first 10% of the chan- 
nel act as reservoir. In the experiment the channel length is 
L x = 444 .4a = 2.0 mm. The strength of the constant driv- 
ing force F ext = Fe x can either be specified directly or by 
definition of the inclination a resulting in F = mg sin a. Un- 
der the influence of external driving the particles still form 
layers. Additionally we observe, both in experiment and in 
simulation, a decrease of the number of layers in the direction 
of motion [25, 66]. The layer transitions are clearly visible 
in Fig. 15, where they are located at x « 420a from 8 to 7 
layers, at x ~ 700a from 7 to 6 layers, and at x ~ 770a from 
6 to 5 layers. 
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FIG. 16: (a) Experiment: Non reworked video microscopy snapshot 
of colloidal particles moving along the lithographically defined chan- 
nel. The channel partition shown has the size (692 x 60 )u.m = 
(153.8 x 13.33)<r, and the interaction strength is V w 72. 
(b) Simulation: Snapshots for a channel with ideal hard walls 
[(573.3 x 45 )pm = (127.4 x 10) a 9 T = 640.5], (c) the same 
as in (b) with the particles at the walls (marked green) kept fixed 
[(573.3 x 45 )\im, V = 5026]. The rectangles mark the region of 
the layer reduction. 

The images of Fig. 16 show in enlargement the part of the 
channel near the region of layer reduction being marked by 
the rectangle. The video microscope snapshot of Fig. 16(a) is 
taken from the experiment [25]. The small white spots at the 
particle centers allow for precise tracking of the particle tra- 
jectories with the video microscope. Similar snapshots we get 
from our BD simulations with either co-moving (Fig. 16(b)) 
or fixed edge particles (Fig. 16(c)). In these two subfigures the 
filled circles represent the particles at their real size relative to 
the channel width. For these highly ordered systems the layer 
transitions take place on the scale of only a few particle diam- 
eters. 



B. Density Gradient along the Channel 

The simulation snapshots above are taken after a time long 
enough for the system to reach a stationary non-equilibrium 
situation. Applying the external driving force to the equili- 
brated channel configuration leads to the build-up of a parti- 
cle density gradient along the channel. This is an effect of 
the chosen boundary conditions at the channel entrance and 
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exit, which leads to a pressure difference between both chan- 
nel ends. After about 10 6 BD time steps this density gradient 
does not change any more, which is the signature of a station- 
ary state. The exact origin of the density gradient is given by 
details of the particle-particle interactions in combination with 
the driving force and will be subject of a separate publication. 




x[c] 

FIG. 17: Simulation: Stationary non-equilibrium density histograms 
along the channel for several values of the slope a (T = 533.75 
and n — 0.4 a -2 ). The vertical line at x — 80a marks the right 
end of the reservoir, i.e., the maximum x-value up to where particles 
are inserted randomly. The inset shows the density gradient in the 
interval x G [150,600]cr as a function of the inclination a being 
obtained from linear fits to the density histograms. Here the line 
connecting the data points serves as a guideline to the eye. 

To study the robustness of the formation of the density gra- 
dient and its connection to the layer reduction in our system 
of gravitationally driven particles, we performed simulations 
for a variety of inclinations a = 0.0° — 10.0° keeping the 
overall particle density fixed atn = 0.4a -2 . The resulting 
stationary non-equilibrium density profiles along the channel 
are shown in Fig. 17. They are calculated from histograms 
of the x-positions of 1000 configurations in stationary non- 
equilibrium. A very significant decrease of the local density 
occurs for x > 700a, which is caused by the open boundary 
at x — 800a. The region < x < 80a acts as reservoir, 
where new particles are inserted at random position whenever 
a particle drops out at the end of the channel. To avoid unnec- 
essary high perturbations due to random particle re-insertion 
in the reservoir the channel is closed at x = 0a by a semiper- 
meable ideal hard wall. Note, that with increasing inclination 
a a depletion layer forms within the reservoir which is the 
result of a greater outflow than input of new particles. 

All density profiles show a nearly linear density gradient in 
the interval x G [150,600]<r, which is maximal for a — 0° 
(cf. Inset of Fig. 17). Even at a = 0° a (osmotic) pressure 
difference between both channel ends exists for the boundary 
conditions used, and a small particle flux is induced. For in- 
clinations a > 1.0° the density gradient becomes almost zero. 
For these inclinations the driving force dominates, and we find 
plug flow of the particles without layer reduction. A decrease 
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of the inclination (driving force) gives rise to an increase of 
the density gradient. Under non-plug flow condition we find 
a self-induced arrangement of the particles to a nearly hexag- 
onal lattice and the occurrence of layer reductions with the 
particles moving across. 



C. Dynamical Properties 

1. Drift Velocity 
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It is also interesting to study the average overall drift ve- 
locity as function of the driving force. The result is shown 
in Fig. 18. For a > 0.5° the particle flow is dominated by 
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FIG. 19: Simulation: Drift velocity histograms of the particles in 
different channel regions. The vertical line marks the expected drift 
velocity for non-interacting particles according to the Drude model 
(cf. equation (9)). 
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FIG. 18: Simulation: Average particle drift velocity in the interval 
x £ [200, 770] a as function of the inclination a or equivalently 
the driving force F = mg sin(a)e x . (L x = 800cr, L y = 10a, 
n — 0Aa~ 2 and T = 533.74). The solid line is the expected 
drift velocity for non-interacting particles due to the external driv- 
ing (Drude model). 

the driving force. This is the regime of plug flow, where the 
particles move with 

(^drift/Drude = y Sm« (9) 

as expected for non-interacting particles. Such a dependency 
was formulated by P. Drude [67] for electrical conduction to 
explain the transport of electrons in metals. For a < 0.5° 
the average drift velocity deviates from the expectation of the 
Drude model. Interestingly, for these inclinations the particles 
move faster than expected. The Drude model is based on a 
friction dependent mobility coefficient, only. For inclinations 
a < 0.2° the diffusion behavior of the particles has to be 
taken into account, too. Therefore, the interplay of the small 
drift and of the diffusion behavior gives rise to a change of the 
mobility in ^-direction. 

Particles moving along the channel get accelerated. This 
becomes obvious from Fig. 19, where histograms of the drift 
velocity together with Gaussian fits in different x-regions of 
the channel are plotted. All hydrodynamic interactions are 
neglected. Generally, the particle velocities v x in ^-direction 
are normally distributed about the average drift velocity. For 



FIG. 20: Experiment: Drift velocity histograms of the particles in 
the full field of view. Shown are the histograms for the bulk and the 
edge particles, respectively. The data points of the bulk particles can 
be fitted well by a single Gaussian, whereas the edge particles need 
to be fitted by a superposition of two Gaussian functions. The results 
of the fits are indicated by the lines. 

the angle a = 0.2° the average drift velocity is (^drift) ~ 
0.081 |xm/ s. In the experiment an inclination of a exp = 0.6° 
was chosen, which results in (thrift ) ~ 0.035 [im/s. The 
velocities of the particles in the experiment are lower as com- 
pared to simulations, possibly due to the influence of hydro- 
dynamic interactions. The comparison between edge particles 
and bulk particles shows the effect of layer changes of the 
particles on the velocities. As mentioned in section IV dislo- 
cations are present along the walls. These dislocations lead to 
an increased number of layer changes for the edge particles. 
During the layer transition the particles move in ^/-direction 
rather than x-direction; therefore we expect to see a super- 
position of 2 velocity distributions in ^-direction, one cen- 
tered around zero for particles changing layers and one cen- 
tered around the velocity of the particles in the edge layer. 
Figure 19 shows Gaussian fits for the bulk particles and the 
edge particles. It is apparent that the velocity distribution of 
the edge particles can be fit by a superposition of 2 Gaussian 
fits, resembling the particles changing lanes (around (Jjn/s) 
and moving straight (0.031 jLtm/s). The different velocities of 
bulk and edge layers are caused by the difference in density 
of the layers. The behavior of bulk particles can be fit with a 
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single Gaussian, because the percentage of particles changing 
layers is much lower than in the edge layer as the transition is 
confined to a small region. 



2. Example Particle Trajectories 
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FIG. 21: Simulation: (a) Example particle trajectories which show 
the dynamical rearrangement of the particles crossing the layer re- 
duction zone from 8 layers to 7 layers. Shown are the trajectories 
for the time interval At = 37.5r s (= 5 • 10 5 BD steps), (b) Corre- 
sponding snapshots of the starting and final configuration. The gen- 
eral color coding (see text) is used. Additionally, all the particles 
which trajectories are shown in (a) have been marked in magenta, 
(c) Histograms of the ^-positions within different x-regions evalu- 
ated for 1.5 • 10 6 BD steps. The peaks of the edge particles are trun- 
cated for clarity reason. The system parameters are identical to those 
of Fig. 15. 

The particles flow across the layer reduction zone 
(cf. Fig. 21), whereas the position of the layer reduction zone 
almost remains unchanged. We show in Fig. 21(a) representa- 
tive particle trajectories for a selection of particles. These are 
marked in magenta in the configuration snapshots (Fig. 21(b)) 
at their beginning and the final location of the trajectories. The 
trajectories clearly show that we do not observe plug flow of 
a crystal, but rather a dynamic behavior of particles moving 
in layers and adapting to the external potential. The parti- 
cles move a distance of about 60a whereas the layer transition 
stays located withing x G [390, 400] a. 

The edge particles are pushed against the ideal hard walls at 
y = Oct and y = L y = 10a by the repulsion of the inner parti- 



cles of the channel. This is the reason for their minimal fluctu- 
ations perpendicular to the flow direction. The corresponding 
fluctuations of the non-edge layers are significantly larger, and 
a small increase of the mobility in ^-direction with increasing 
wall separation is found. In the central region some particles 
change very abruptly from one layer to another whereas oth- 
ers shift more smoothly. The particles in the layers next to 
the edge layers only show a small and smooth change in their 
^/-position. In the regions with fixed number of layers no parti- 
cle transitions between layers are observed for our simulation 
parameters. 

All particles are identical. In Fig. 21(b) we just color coded 
the particles according to the number of nearest neighbor par- 
ticles they have. Bulk particles with six nearest neighbors 
and all edge particles are marked blue, whereas red particles 
have a fivefold symmetry and green particles have a sevenfold 
symmetry of nearest neighbors. The actual number of nearest 
neighbors is determined using a Delaunay triangulation. In 
the start configuration three defect pairs (dislocations) are in 
the region of the layer transition form 8 to 7 layers, whereas 
in the final configuration the layer reduction position is con- 
nected to a single dislocation. The slightly higher density of 
the edge particles gives rise to the scattered green particles in 
the next edge layer. 

For the same system we analyze the density profiles trans- 
verse to the walls within several sub-regions along the chan- 
nel. Therefore we evaluate 1.5 • 10 6 BD steps correspond- 
ing to a time interval of At w 122. 5r^. The full density 
profile for x G [100, 740] a (black curve) is a superposition 
of several profiles connected to distinct layering. Highly or- 
dered layer structures with sharply peaked density profiles oc- 
cur for eight layers in x G [100, 380] a (red curve), seven 
layers in x G [440, 630] a (blue curve), and six layers in 
x G [670, 740] a (cyan curve). The x-regions in between are 
the layer transition regions. 

In Fig. 22 we explicitly plot the superposition of 191 1 video 
microscopy snapshots of the experimental system and 3000 
configuration snapshots used for the density profile evalua- 
tion above. In the experiment the particles move on average 
(Ax) w 670a. The layer reduction zone is confined in the 
interval x G (5, 50) a. In the simulation the particles have 
moved forward on average the distance (Ax) « 202a, i.e. 
more than a quarter of the channel length. The layer transi- 
tion positions remain located within an interval of length 45a. 
The particles are inserted at a random position in the region 
x G (0, 80) a. Perturbations of the configuration due to the 
random particle insertion heal after a few BD steps. There- 
fore the configuration for x > 90a is not influenced by this 
particle re-insertion method. 



3. Defect Removal 

Sometimes "defects" remain after the point of layer reduc- 
tion, which vanish on further flow. Here we call a defect a 
pair of particles having 7 and 5 neighbors respectively which 
disturb a given layer configuration. These can be identified 
from dips they form in the local layer order parameter defined 
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4. Diffusion Behavior 




FIG. 22: Superimposed configuration snapshots: (Top) of the ex- 
periment, (Bottom) of the simulation for 1.5 • 10 6 BD steps, which 
corresponds to At « 122. 5r s . (B = 0.5 mT, T = 533.74, 
L x = 800a, L y = 10a, n = 0Aa~ 2 , a = 0.04°) 
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FIG. 24: Particle diffusion behavior in the experiment. 

In section IV the diffusion behavior of particles in the equi- 
librated channel was discussed. The question arises whether a 
similar behavior can be observed for the driven system. Fig- 
ure 24 shows the MSD in a> and ^/-direction for one of the bulk 
layers in the experiment after the driven motion of the system 
has been subtracted. In x-direction we see a linear increase 
of the MSD at short times, which starts to saturate at longer 
times. A clear transition to SFD, as it was found in the simu- 
lation data for the equilibrated channel, cannot be found. One 
reason can be given by the rather short times at which the ex- 
perimental MSD could be obtained. In ^-direction we see sat- 
uration at rather short time scales. This behavior is due to the 
fact that the particles move in stable layers during most of the 
experiment and are therefore restricted in their ^-movement. 
A similar behavior was shown in Fig. 14 for a channel width 
which induces stable layers. 



FIG. 23: Simulation: The sequence of configuration snapshots (a) - 
(h) shows the process of a vanishing "defect" after the layer transition 
zone (marked by the black rectangle) due to the change of particle 
marked in orange into the edge layer. The snapshots have been taken 
every 500 BD time-steps, i.e. At = 0.0375r s . 



in equation (5) of the current configuration. Generally, small 
density gradients along the channel give rise to a larger num- 
ber of defects than higher density gradients. This already is a 
hint on the close connection of the occurrence of layer transi- 
tions to the local number density. A defect can be neutralized 
by a particle changing into the edge layer. Such a neutraliza- 
tion process of two defects is shown in the sequence of config- 
uration snapshots of Fig. 23 taken every 500 BD time-steps. 
The orange colored particle moves into the edge layer and 
thereby removes the perturbation of the layered structure of 7 
layers after the position of the layer transition region marked 
by the black rectangle. In the final snapshot 23(h) seven un- 
perturbed layers remain. Recognize that again the x-position 
of the layer reduction remains unchanged. 



D. Connection between the Layer Transition and the Density 
Gradient 

The reduction of the number of layers originates from a 
density gradient along the channel. The local particle density 
p{x) inside the channel is shown in Fig. 25 together with the 
local lattice constants d x and d y . The particle separations of 
neighboring particles in x- and ^/-direction are used to calcu- 
late the local lattice constant d of the triangular lattice. Due to 
the density gradient along the channel, the ordered structure 
is not in its equilibrium configuration at all points along the 
channel. Thus the local lattice constant d x , calculated from 
the particle separations in ^-direction, can deviate from the lo- 
cal lattice constant d y , calculated from the particle separations 
in ^-direction and multiplication with the factor 2 / >/3. At the 
left end of the channel, d x increases to larger values than d y , 
indicating that the ordered structure is stretched along the x- 
axis. At the position of the layer reduction the system changes 
back to a situation, where d x is smaller than d y by decreasing 
d x and increasing d y by about 20% simultaneously. These 
changes of separations compensate each other and result in a 
continuous change in the local density at the position of the 
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Static Stretching Analysis 



layers 
7 layers 
6 layers 
5 layers 



FIG. 25: Simulation: (a) Local lattice constants d x and d y and local 
particle density p, (b) Corresponding local layer order parameters 
flayer, n L ■ The system parameters are: L x = 800a, L y = 10a, 
n = o'.4cr- 2 , and r = 533.74. 



layer reduction. The behavior of the system shows that the 
stretching of the ordered structure before the layer reduction 
causes an instability towards decreasing the number of layers. 
This decrease compresses the system along the x-direction, 
but apparently lowers the total energy of the system. 
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FIG. 27: Result of the stretching analysis of static channel configu- 
rations of a channel with the width L y = 10 and dipolar pair interac- 
tion: Shown are the potential energies per particle of different layer 
configurations as function of the particle density p. 

The scenario above can be qualitatively confirmed by the 
following rough estimation: Starting from an ideal triangular 
configuration with a given number of layers (n/) in a channel 
of fixed width, we calculate the potential energy per particle 
for different particle densities by scaling the channel length 
of the static configuration only. Plots of these energies per 
particle for different values of n\ as function of the particle 
density p are shown in Fig. 27. The intersection points, which 
are determined from linear approximation of both curves in 
the region of intersection, serve as a rough estimate of the 
density at the layer transition point. For the given system the 
values for the transition 8^7 layers are: ps^7 ~ 0.467a -2 , 
and for the transition 7^6 layers: ^7^6 ~ 0.345a -2 . The 
full circles mark the perfect triangular lattices configurations 
of the respective number of layers. So, we can conclude that a 
given layer structure is stable for up to slightly overstretched 
perfect triangular configurations. 

They show clear intersection points, indicating that for a 
stretched configuration with m layers in x-direction it can be- 
come energetically more favorable to switch to a compressed 
configuration with (ni — 1) layers. 



Equilibrium Configurations for Confinement with Non-Parallel 

Walls 



FIG. 26: Simulation: Local lattice constant d y as function of the 
local particle density p for various inclinations a. 

Layer transitions occur at almost identical values of the lo- 
cal particle density for various inclinations as can be seen in 
Fig. 26. Here the local lattice constant d y (x) is plotted as a 
function of the local particle density p(x). Transitions from 8 
to 7 layers occur when p(x) becomes smaller than 0.42, tran- 
sitions from 7^6 layers for p(x) < 0.3, and transitions 
from 6^5 layers for p(x) < 0.21. 



Also equilibrium BD simulations, i.e. simulations with no 
external driving force (Ff xt = 0), of closed channels with 
non-parallel walls result in a density gradient in the direc- 
tion of decreasing channel width. Here, confinement induced 
arrangement of the particles into different number of layers 
takes place. The particles just fluctuate about their equilib- 
rium positions. A snapshot of such an equilibrium config- 
uration is shown in Fig. 28 where the confining funnel has 
the small opening angle apunnei = 0.143°, i.e. over the full 
channel length of L x = 800a the channel width decreases 
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FIG. 28: Simulation: Equilibrated configuration snapshot of a fun- 
nel geometry with opening angle chopper = 0.143°. All walls are 
modeled as hard walls. No driving force is applied and the particles 
interaction is dipolar (T = 625.125). 



by AL y = 2(j. This kind of layer transition is a purely ge- 
ometrical effect, whereas in the case of parallel walls and a 
constant longitudinal driving field the occurrence of the den- 
sity gradient has a dynamical origin. In both cases the number 
of layers which form depends on the value of the local particle 
density p(x). 
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FIG. 30: Snapshots of defect configurations obtained from a Delau- 
nay triangulation of the particles moving in the channel. The parti- 
cles are coded according to the number of their nearest neighbors. 
Open circles mark the bulk particles with 6 nearest neighbors and 
the edge particles, symbol x corresponds to a fivefold symmetry, 
and symbol V to a sevenfold symmetry, (a) Experiment: In order to 
minimize the effects of fluctuations on a short time scale, 50 images 
have been averaged, (b) BD simulation for a channel with parallel 
walls. 



E. Comparison with the Experiment 

The experimental result of the density gradient as well as 
the interparticle distances are shown in Fig. 29. The behav- 
ior closely resembles the behavior of the simulated system 
(cf. Fig. 25(a)). The distance in x-direction, d x , is continu- 
ously stretched while the distance in ^-direction increases in 
a sharp step at the position of the layer reduction. The density 
decreases monotonously along the direction of motion of the 
particles by about 20%. 
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FIG. 29: Experiment: Local lattice constants d x and d y and local 
particle density. The results are obtained for the systems of which is 
shown in Fig. 16(a). 

In Fig. 30 snapshots of non-equilibrium defect configura- 
tions are shown both for the experiment (a) and for the sim- 
ulation (b). They reveal that the system is nearly triangular 
left and right of the point of layer reduction. The change is 
marked by a single defect only. The number of layers is re- 
duced one by one. Reductions of two or more layers have 
not been observed in experiment or in simulation. Naturally, 
this reduction produces a defect at the point of the transition. 
Since the position of the layer reduction is mainly determined 



by the density gradient, its location remains stable with time 
on average. A more detailed analysis reveals, however, that 
the transition point oscillates back and forth around this aver- 
age position. At the transition the driven particles in the bulk 
layers have to change the layer, causing the transition to move 
a little bit in direction of the flow. A particle changing into the 
edge layer can neutralize the defect of the transition locally. 
This causes a reconfiguration of the ordered structure, which 
in turn gives rise to repositioning of the layer reduction zone 
back to a region of higher density. 



F. Oscillatory Behavior of the Layer Transition 
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FIG. 31: Simulation: Movement of the x-position of layer transition 
for the transitions 8^7 layers and 7^6 layers. The system 
parameters are identical with those of Fig. 22. 

There are various ways of numerically localizing the po- 
sition of the layer transition. One can either make use of 
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the clear discontinuity of the local layer order parameters 
^ layer, m(x,t) (cf. equation 5) appropriate for the transition 
from ni to n\ — 1 layers, or of the location of the discontinuity 
of the local lattice constant d y {x,t). The local orientational 
order parameter which is often used for 2D systems [18], 
is not so significant for this system, as it is very sensitive 
to any perturbation of the sixfold symmetry. The first three 
methods have been used to study the position of the transition 
from 8 to 7 layers. The result is given in Fig. 31. A com- 
parison of all three methods mentioned is given in Fig. 32 for 
the transition from 8 to 7 layers. The local order parameters 
are calculated within bins of size l x = 2a in flow direction 
limiting the x-resolution. As can be seen all methods give 
similar results, but special care needs to be taken in the pres- 
ence of defects, which occur close to the layer transition for 
t > 120r B = 25.2 • 10 3 s. 
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FIG. 32: Simulation: Movement of the x-position of layer transition 
for the transitions 8^7 for an inclination a — 0.2° and otherwise 
identical parameters as in the previous figure. 



In the non-equilibrium steady state situation the position of 
the layer reduction zone oscillates about a certain x-position. 
This can also be seen in the experimental data, as shown in 
Fig. 33 (evaluated from the discontinuity in d y ). 
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FIG. 33: Experiment: Movement of the x -position of layer transition 
for the transitions 8^7. 



G. Influence of the Particle Interaction Range 

In order to study the influence of the particle interaction 
range, we implemented the screened Coulomb (YHC) pair in- 



Vij(rij) = < 



oo 
I 



exp(-K D (nj - a)) 



nj < a 

a < nj < r cut 



(10) 

with the inverse Debye screening length kjj which interpo- 
lates the potential between the hard core case (for kd — » oo) 
and the unscreened Coulomb potential (for k,d =0). Vo is the 
value of the pair potential at contact which can be written as 



PV = 



X? 



(i + n D (j/2y 



where Z is the charge of the colloids and \b = 
e 2 /(4:7reoe s kBT) is the so-called Bjerrum length of the sol- 
vent with permittivity e s . 
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FIG. 34: Simulation: (a) Local lattice constants d x and d y and local 
particle density p in the BD simulation of a system with screened 
Coulomb interaction, (b) Corresponding local layer order parameters 
flayer, n t ■ The system parameters are: L x = 800cr, L y = 10cr, 
n = 6.4cr" 2 , pV = 400, k, d = 4.0a" 1 , Tyhc = 448.4, and 
a = 0.2°. 

Figure 34 is the analogous plot to Fig. 27 for a system 
of YHC-particles with the contact value (3Vq — 400 and 
hv£) = 4. Oct -1 . Under these simulation conditions no layer- 
transition as for the dipolar system is found. For x > 450a 
the particles are ordered in 7 layers, but for smaller values 
only a few islands of particles arranged in layers can be iden- 
tified from the local order parameters along the channel in 
Fig. 34(b). The interaction range of a YHC system with 
kd = 4. Oct -1 is much smaller than for the dipolar system, 
because of the stronger decay of the pair potential. This decay 
is also the reason for the large fluctuations of the local lattice 
constant d x in Fig. 34(a). A density gradient can not form 
along the channel, and so no layer transition is found. The 
particles need to be strongly coupled with their neighboring 
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particles to form a density gradient, i.e. the pair interaction 
range has to be at minimum of the order of the average parti- 
cle spacing. 



VI. ALTERNATIVE BOUNDARY CONDITIONS IN FLOW 
DIRECTION 

The connection of the channel to the two reservoirs has 
great influence on the characteristics of the stationary non- 
equilibrium density profile along the channel. Therefore, we 
performed simulations with an alternative boundary condi- 
tion, where the constant external driving force only acts within 
the interval x G [100, 700] a and a periodic boundary condi- 
tion is applied in x-direction. Figure 35 shows the resulting 
stationary non-equilibrium density profiles after 3 • 10 6 BD 
time steps for a selection of inclinations a of a dipolar system. 
For each inclination two curves are plotted which correspond 



to the two channel widths L lt 



8a and L ? , 



10a. Obvi- 



ously, the steady- state density profile along the channel does 
not depend on the channel width. 
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In the stationary non-equilibrium state the density profile in 
the reservoirs can be approximated by a linear gradient. The 
net flux J in the reservoirs fulfills Fick's law 



J 



k B T 
2l 



(pi ~ Po) 



(11) 



where po = p( x — 100a) and p\ = p(x = 700a) are 
the local number densities at the channel beginning and end 
respectively. Due to the periodic boundary condition in in- 
direction this is equal to the net flux in the channel region 
x G [100, 700)a. Therefore, J may be approximated by the 
slope of the linear density profiles in the two reservoir regions. 
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FIG. 35: Simulation: Density profiles for a selection of inclinations 
a of a system with the inclination (i.e. the driving force) applied 
only within the region x G [100, 700] a. The system is periodic in re- 
direction. Shown are histograms obtained by evaluation of 1000 con- 
figurations of the system having reached a stationary non-equilibrium 
situation (after w 2 • 10 6 BD time steps). The applied magnetic 
field strength is B = 0.25 mT, V = 133.4, and the overall parti- 
cle density n — 0.4a -2 . For better clarity, we replicate the interval 
x G [0, 100] a again on the right hand side of the diagram. 

Comparison of these density profiles with those of Fig. 17 
highlights the strong influence of the different realization of 
the reservoirs. All simulations are started from a homoge- 
neous particle distribution of local density p = 0.4a -2 . In- 
stead of a density decrease we find in Fig. 35 a buildup of 
the local density occurring due to the filling of the reservoir 
at the channel end. This corresponds to the experimental sit- 
uation, where the reservoir at the channel end is filled. For 
the small inclination a = 0.04° a linearly increasing den- 
sity profile is obtained within the channel region. Higher in- 
clinations lead to deviation from such a linear profile. For 
a = 0.2° a constant profile with local density p « 0.275a -2 
in x G [100, 400] a is followed by a sharp increase of the local 
density up to p = 0.67a -2 at the channel end at x = 700a. 



FIG. 36: Simulation: Stationary non-equilibrium situations of sys- 
tems where the driving force is applied only in x G [0, 700] a 
for a selection of inclinations: (a) a = 0.04°, (b) a = 0.08°, 
(c) a = 0.1°, (d) a = 0.2°. For every inclination we show the av- 
erage local layer order parameters ^i ay er,n z and the corresponding 
superposition of 1000 snapshots. The other simulation parameters 
are: L x 800a, L y 10a, n 0.4a -2 , B 0.25 mT, and 
T 133.4. 

Figure 36 shows the layer order parameters ^i ay er, m for 
a selection of inclinations a in combination with the corre- 
sponding superimposed configurations. Clearly, the layer con- 
figuration and the number of layer transitions can be tuned 
by the strength of the driving force for the realization of the 
boundary condition 2 of the flow. Increasing a leads to multi- 
ple transitions. Interestingly, the layer transitions from 7 to 8 
layers occur at identical ^-positions in the figures 36(b)-(d). 
As before, the particle flow across the position of the layer 
transition, which remains fixed in position. 

Systems with Screened Coulomb Interaction 

In Fig. 37 we show the equilibrium density profiles trans- 
verse to the confining walls of a YHC system for a selection 
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of values and the value at particle contact /3Vq = 50. The 
total particle density is n = 0.45a -2 which corresponds to a 
packing fraction j] w 0.79. Forft£> = 2a -1 boundary induced 
layering is found which becomes less pronounced for increas- 
ing k, i.e. decreasing interaction range. For > 4a -1 the 
systems are fluid in the equilibrium state at this packing frac- 
tion, and only a depletion layer between the edge and the bulk 
particles can be seen. 




niinet 1 ^] 



urations with a time separation of At = 500 BD steps af- 
ter 1.4 • 10 6 BD steps for the case of the alternative boundary 
condition in flow direction. The driving force corresponding 
to an inclination of a = 0.1° acts within x G [100, 700] a. 
All four superimposed configurations show the formation of 
layers near the channel end at x = 700a, where the particles 
enter the reservoir. In Fig. 38(a) the characteristic interaction 
range of the YHC pair-potential is greater than the average 
particle spacing R. For this case we find multiple layer tran- 
sitions from 5 layers up to 8 layers along the channel. The 
system behavior is similar to the situation of the dipolar sys- 
tems. With increasing values of k,d less layer transitions are 
observed. Figures 38(b)-(d) show increasing depletion zones 
at the channel start at x = 100a. These depletion zones are 
followed by regions where the particles are in the liquid state. 
Notice, that for Kb = 8a -1 and kd = 12a -1 the systems are 
in the liquid state in equilibrium, too (cf. Fig. 38). The corre- 
sponding density profiles in x-direction are given in Fig. 38. 



FIG. 37: Simulation: Full density profiles transverse to the confining 
walls for L y — 8a of systems with screened Coulomb interaction 
for a selection of interaction ranges kd- The contact value of the 
potential is /3Vo = 50. 

The average particle separation of the unbounded system 
has the value R w 1.38a. For > 4a -1 the characteristic 
interaction range is a + k^ 1 < 1.25a which is smaller than 
R. 
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FIG. 38: Simulation: Superimposed configurations of systems with 
screened Coulomb pair interaction for a selection of inverse screen- 
ing lengths: (a) kd = 2a~ 1 , (b) kd = ^o~ A ", (c) kd = 8<j _1 , 
and (d) «d = 12<j _1 . The particle transport is induced for x G 
[0, 700] a by the inclination a = 0.1°. 

Now, we plot in Fig. 38 the superposition of 100 config- 
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FIG. 39: Simulation: Density profiles along the channel for a se- 
lection of Debye screening lengths kd of a YHC system ((3Vo = 
50). The driving force is applied only within the channel region 
x G [100, 700] a and the system is periodic in x-direction. These 
profiles correspond to the superimposed configurations of Fig. 38. 

The systems with kd — 8a _1 and Kb = 12a -1 show a 
rapid increase of the local density from about 0.5a -2 up to 
values greater than 0.8a -2 in the interval x G [600, 700] a. 
The particles under the influence of the constant driving force 
are blocked due to filling of the reservoir at the channel end. 
During the simulation run the particles pile up at the interface 
between the channel and the reservoir, because the particles of 
the channel are pushed into the reservoir but within the reser- 
voir the particles diffuse almost freely due to the short range 
of the YHC interaction (high values of kd)- This leads to 
a situation where the influx into the reservoir is greater than 
the particle drift within the reservoir being the reason for the 
sharp density gradients, which lead to the sudden onset of a 
layered structure with 8 layers in the figures 38(c)-(d). For 
KjB = 12a -1 even a layer transition to 9 layers takes place 
due to local density values greater than 0.9a -2 which is not 
observed for the other three cases. Alternatively, the particle 
flux can be blocked in a controlled fashion by creating so- 
called laser barriers perpendicular to the driving field, as we 
will show in the following section. 
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VII. CONCLUSION 

We have reported on a variety of ordering and transport 
phenomena which are induced by the confinement of colloidal 
particles to microchannels and by the application of a constant 
driving force along the channel. We have analyzed the particle 
behavior both under equilibrium and under (stationary) non- 
equilibrium conditions both in experiment and by Brownian 
Dynamics simulations. 

First, we have studied the self-assembly of repulsive parti- 
cles under equilibrium conditions, i.e. without a driving force 
applied.We have observed a boundary induced formation of a 
global layered structure in the channels. Such a behavior is 
known for a variety of related systems [14, 22, 31, 34, 42]. 
Systematically, we have analyzed the influence of the chan- 
nel width L y and the influence of the strength of the dipolar 
particle repulsion. Based on the order parameter we have cal- 
culated the phase diagram of laterally confined superparamag- 
netic particles as a function of the channel width L y within the 
solid state. We have observed a re-entrant behavior as a func- 
tion of L y , where the system behavior oscillates from solid- 
like to liquid-like. When the channel width is increased, a 
periodic destabilization of the layered structure with n\ layers 
takes place, and the system switches to a structure with n\ + 1 
layers. The bulk defect concentration C^ efect shows periodic 
oscillations as a function of the channel width L y , but not as 
a function of the dimensionless interaction parameter T. The 
period of the oscillations is ~ R, where R denotes the average 
distance of two neighboring particle layers in an unbounded 
hexagonal system. Such a behavior previously was reported 
as a result of both experiments and simulations of a similar 
system by Haghgooie et al [22, 23, 30]. Our data show excel- 
lent qualitative and quantitative agreement with their results. 

For very small channel widths L y < la, where the parti- 
cles cannot pass each other, we have observed the well known 
single file diffusion behavior [44, 45, 47]. Our model sys- 
tem allows for a systematic analysis of the diffusion behavior 
of layered structures. In first studies, we have compared the 
longitudinal and transversal particle diffusion behavior for the 
channel width L y = 10a, where the system globally forms 
7 layers, to a channel with L y = 9a, where no globally 
layered structure exists in the solid state. The diffusion be- 
havior at intermediate time scales is very different for both 
cases. In the presence of global layers, the transversal mean 
square displacement (Ay 2 ) has a constant plateau. Addition- 
ally, at longer times a deviation from the Fick diffusion be- 
havior (anomalous diffusion) is expected by our first simula- 
tion results. Evaluation of the experimental data confirm this 
deviation from the Fick diffusion behavior. These effects are 
absent for the system with L y = 9a. Therefore, it will be in- 
teresting, to do a systematic analysis of this diffusion behavior 
in the future. 

We have predicted and systematically analyzed the phe- 
nomenon of particle layer reduction under the influence of 
a constant driving force acting along the channel. For small 
driving forces F ext , where the particles are not yet in the 
regime of plug flow the superparamagnetic particles dynami- 
cally re-arrange into different numbers of layers during trans- 



port through the channels. We have found, that along the 
channel the number of layers decreases gradually by steps of 
one. The occurrence of the layer reduction has been confirmed 
by the experiments. In the experiments, the massive particles 
sediment to the bottom of the channel due to gravity, and there 
they form a quasi-2D system. After having equilibrated the 
system, the whole setup is tilted, so that the colloidal particles 
are driven through a lithographically fabricated microchannel 
under the influence of gravity. 

In very good qualitative agreement with the experiments we 
have shown that the reduction of layers originates from a den- 
sity gradient along the channel. Quantitative differences are 
expected, because the Stokes diffusion coefficient Do, which 
is valid for unbounded systems and is used in the simulations, 
differs from the real diffusion coefficient in presence of the 
confinement of the experimental setup [68]. 

The reduction of layers takes place for specific values of 
the local density p(x) and within a distance of only a few par- 
ticle diameters. We have explicitly shown that the particles 
flow across the regions of layer reduction and thereby dynam- 
ically adjust to the local density p(x). The origin of the local 
density gradient is not fully understood yet. But additional 
simulation studies of systems with screened Coulomb particle 
interaction, where the interaction range has been varied, have 
shown that a longitudinal density gradient and consequently 
layer transitions occur for particle interaction ranges which 
are greater than the average distance of the particles from their 
neighbors. For particle pair potentials with smaller interac- 
tion ranges than the average nearest neighbor separation we 
observe that the layer transition region smears out, because 
more particle defects occur due to a smaller density gradient. 
No layer transitions will be observed for the model-case of 
hard-core particles. For our choice of boundary conditions 
we have found, that the density gradient becomes more pro- 
nounced with decreasing inclination a, i.e. with decreasing 
driving force. The density decrease is maximum at a = 0.0°, 
because the particle re-insertion scheme, which we used, in- 
duces a pressure difference between both channel ends, even 
in the case when no external driving force has been applied. 

Generally, we have seen both in simulations and in experi- 
ments that the local density decreases monotonically and con- 
tinuously along the channel. In front of a layer transition the 
local structure is stretched in longitudinal direction, whereas 
after the layer transition the structure is longitudinally com- 
pressed and one layer has disappeared. Therefore, the local 
lattice constant d x {x) in longitudinal direction increases up to 
the position of the layer transition, at which it shows a non- 
continuous decrease. Simultaneously, the local lattice con- 
stant d y (x) in transversal direction is constant in front of the 
layer transition, at which it jumps to the next level according 
to the number of layers and remains constant again. Both ef- 
fects compensate each other and thus explain the continuous 
behavior of the local density along the channel. 

By a static stretching analysis we have confirmed that a cer- 
tain layered structure becomes energetically unstable and thus 
changes to a structure, where it has one layer less. The esti- 
mated values of the local density, where the transition takes 
place, are in quite good agreement with the observation. In 
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stationary non-equilibrium the position of the layer transition 
oscillates about a fixed position. The amplitude of the oscil- 
lations depends on the strength of the particle interaction. We 
have shown, that the oscillations of the layer transition can 
either be analyzed by the appropriate local layer order param- 
eters ^iayer,n z (%) or by the local lattice constant d y (x). 

Each layer transition is connected to a defect, which is de- 
fined by a pair of particles with five and seven nearest neigh- 
bors respectively. Additional periodic defects have been ob- 
served along the channel walls. Due to the purely repulsive 
particle interaction the edge particles are pushed against the 
flat walls. This leads to very small transverse fluctuations of 
the edge particles and a slightly higher line density of the edge 
particles than of the particles belonging to the layers in the 
central region of the channel. 



It has been shown, that channel walls made of periodically 
fixed particles give rise to shear effects between the particles 
of the central layers, which move faster, than the particles, 
which are in the layer next to the edge particles. The latter 
particles show small oscillations about the average drift ve- 
locity. 

The results shown concern a rather simple classical model 
system. The observed phenomena, however, will take place 
in any systems of long range interacting particles which are 
driven through a constriction. Therefore the results which 
have been gained from the studies of this system can be seen 
as a first step in the understanding of transport processes in 
many biological and quantum systems. 
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SFB TR6 and the NIC, HLRS, and SSC. 
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